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SUMMARY 

This research has addressed the problem of real time guidance and optimal control of 
Aeroassisted Orbit Transfer Vehicles (AOTV's), using singular perturbation theory as an 
underlying method of analysis. Trajectories have been optimized with the objective of 
minimum energy expenditure in the atmospheric phase of the maneuver. Two major problem 
areas were addressed: optimal reentry, and synergetic plane change with aeroglide. For the 
reentry problem, several reduced order models were analyzed with the objective of optimal 
changes in heading with minimum energy loss. It has been demonstrated that a further model 
order reduction to a single state model is possible through the application of singular 
perturbation theory . The optimal solution for the reduced problem defines an optimal 
altitude profile dependent on the current energy level of the vehicle. A separate boundary 
layer analysis is used to account for altitude and flight path angle dynamics, and to obtain lift 
and bank angle control solutions. By considering alternative approximations to solve the 
boundary layer problem, three guidance laws were derived, each having an analytic feedback 
form. The guidance laws were evaluated using a Maneuvering Reentry Research Vehicle 
model and all three were found to be near optimal. 

For the problem of synergetic plane change with aeroglide, a difficult terminal 
boundary layer control problem arrises which to date has been found to be analytically 
intractable. Thus a predictive/corrective solution was developed to satisfy the terminal 
constraints on altitude and flight path angle. A composite guidance solution was obtained by 
combining the optimal reentry solution with the predictive/corrective guidance method. 
Numerical comparisons with the corresponding optimal trajectory solutions show that the 
resulting performance is very close to optimal. 

An attempt was made to obtain numerically optimized trajectories for the case where 
heating rate is constrained. A first order state variable inequality constraint was imposed on 
the full order AOTV point mass equations of motion, using a simple aerodynamic heating rate 
model. For high heating rate limits (just below the peak heating rate for the unconstrained 
case), the resulting solution appears to satisfy the first order necessary conditions for a 
"touch point" problem, where the constraint is met at a single point. Lower heating rate 
limits likely result in’ a constrained arc, of finite duration. Unfortunately, numerically 
converged optimal trajectories for this range of solutions could not be obtained using the 
multiple shooting method employed in this research. 
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SECTION 1 
INTRODUCTION 

Energy state approximations combined with singular perturbation theory 
have proven useful in aircraft trajectory optimization, both in obtaining 
algebraic control solutions and in satisfying trajectory and control 
constraints [1-4]. However, the underlying flat earth and constant 
gravitational field assumptions in aircraft modeling do not apply to 
hypersonic vehicles. Moreover, the use of singular perturbation theory 
requires an inherent time scale separation in the problem formulation for 
successful application. The intent of this research effort has been to 
explore the usefulness of singular perturbation analysis in the development 
of real time guidance algorithms for problems related to AOTV maneuvers. 

The problem of optimal atmospheric heading change with minimum energy 
loss has application to maneuvering reentry vehicle guidance and to 
aeroassisted orbit transfer vehicle (AOTV) guidance. The problem of 
aeroassisted orbit plane change requires the use of three impulses - one to 
deorbit, one to reorbit and one to recircularize at the new orbit. The orbit 
plane change is effected entirely in the atmosphere through the use of lift 
and bank angle control. Circular orbit plane changes in which the initial 
and final orbital altitudes are equal were studied in [5-9]. These studies 
considered various problem formulations with the underlying approximation 
that an expression related to the sum of the centrifugal and gravitational 
forces (Loh's term) is constant or piecewise constant over the atmospheric 
maneuver. Furthermore, in the absence of heating constraints, the optimal 
trajectories are of short duration, and the inclination change is closely 
approximated by the heading change. For this situation, the dynamics can be 
reduced to fourth order, and minimization of fuel consumption is closely 
approximated by minimizing the energy loss in the atmospheric portion of the 
trajectory [6]. As a point of reference, the optimal AOTV maneuver requires 
approximately 50% of the fuel needed for the single impulse pure propulsive 
maneuver in the case of a 40 degree low Earth orbit plane change. In [10], a 
regular perturbation method is used to remove the approximations related to 
Loh's term in the earlier work and demonstrates a significant improvement 
over the solutions in [5-9]. However, this approach requires a quadrature at 



each update of the control solution and the approach can not be readily 
extended to- include the effect of heating constraints. 

Examples of numerical optimization studies related to orbit plane change 
can be found in [11-13]. In particular, [12-13] examine the effect of a 
heating rate constraint with thrusting in the atmosphere. Since the duration 
of the AOTV maneuver is much greater when a heating rate constraint is 
enforced, it is necessary to consider a more complete set of dynamics which 
includes the cross range angle. Ref. [12] treats the problem of optimal 
aerocrui’se (flight at constant altitude and velocity), and does not consider 
the transitions to and from the cruise condition. In [13], a more general 
problem is treated with a constant thrust segment inserted during the 
atmospheric phase. Thus, in our work we decided to place emphasis on optimal 
aeroglide (no thrusting in the atmosphere) subject to a heating rate 
constraint. A more complete account of related work in noncoplanar transfer, 
including other competing transfer modes, can be found in [14]. 

Section 2 of this report presents the problem formulation and issues 
related to model order reduction. Section 3 treats the problem of optimal 
heading change with minimum energy loss in the context of singular 
perturbation analysis. Section 4 addresses the AOTV synergetic plane change 
problem by introducing a predictive/corrective solution to satisfy the 
terminal constraints on altitude and flight path angle. Section 5 summarizes 
the results for a numerical study of the effects of a heating rate constraint 
on the AOTV synergetic plane change problem. Section 6 summarizes the 
results and recommendations for further research along this line. 



SECTION 2 

PROBLEM FORMULATION 


2.1 Equations of Motion 

The equations of motion for gliding flight about a sperical nonrotating 
Earth are given by: 


d9/dt = Vcosycosi|)/rcos<j): r = r s +h 

(1) 

dc^/dt = Vcosysin^/r 

(2) 

dh/dt = Vsiny 

(3) 

2 

dV/dt = -D/m -gsiny: g = v/r 

(4) 

2 

Vdy/dt = Lcosp/m + (V /r-g)cosy 

(5) 

2 

V dip/dt = Lsinp/mcosy - (V /r )cosycos<|>tan<j> 

(6) 


where 0 is the longitude, <{> is the latitude, h is the altitude, r $ is the 

Earth's radius, V is the velocity, y is the flight path angle, f is the 

heading angle, D is the drag force, v is the Earth's gravitational constant, 
and m is the vehicle mass. The control variables are the lift force (L) and 
bank angl e (p) . 

The orbit inclination angle is given by the relation 

cos i = cos<j>cos<j> (7) 

The plane change is the angle between the normals to the initial and 

final orbital planes. The actual inclination (i) is defined relative to the 

equatorial plane. Many studies on optimal plane change have taken the 

equatorial plane as the initial plane, in which case the plane change equals 
the final inclination angle (i^). However, it has been shown in [12] that, 
under the assumption of spherical symmetry, maximizing the inclination is 

equivalent to maximizing the plane change angle provided that the deorbit 

burn is properly timed so that the plane change occurs at the proper 

location. This depends only on the location of the ascending node for the 
final orbit plane. 



Thus, there is no loss in generality in assuming that the initial plane 
is the equitorial plane, and that the inclination change is the final 
inclination angle. A second consequence of th*is fact is that 0 becomes an 

ignorable coordinate in most AOTV optimization problems since it does not 
appear in the right hand side of the equations of motion, and it does not 
enter the boundary conditions or the performance index for optimal control 
problem formulations of practical interest. 

For short duration maneuvers, the cross range angle can be treated as 
being negligibly small in ( 1 ) , ( 6 ) and (7). In this case the inclination 
change is approximated by the heading change, and <j> also becomes an ignorable 
coordinate. Thus it is possible to reduce the equations of motion to a four 
state model, which for the purposes of this study are expressed in the 
following form: 


V 

dh/dt = Vsiny 

(8) 


dE/dt = -C*S(l+A 2 )pV 3 /4m 

(9) 

w 

dy/dt = (C*pSV/2m) (Acosy+Mcosy) 

(10) 


* 

dif»/dt = C. ApSVsi ny/2mcosy 

(ID 


where 


C D = C Do + KC L 


( 12 ) 


M(h,V) = (2m/C L $)[l-v/rr]/ P r, r = r g +h (13) 

In these equations the superscript * denotes the maximum lift-to drag values: 

il/2 r * = ?r (14) 


C L = (C D0 /K)- 


C D = 2C Do 


and A is the normalized lift coefficient 


A = C L /C L 


(15) 
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Here we have employed E as a state variable in place of velocity (V), where 




E = V 2 /2 - v/r 


(16) 


In [6] velocity is used as a state variable, and the gravity component was 
ignored in the velocity rate equation. One advantage to using E as a state 
variable is that (9) is independent of y. In (16), the reference point for 
zero potential energy is taken at r = ». This transformation Implies that 
wherever V appears in the equations, it is replaced by [2(E+v/r)] ' . The 
control variables are X and the bank angle (y). Under the hypothesis that 
the cross range angle is small, closely approximates the. change in orbit 
inclination. 

A further reduction to a third order model is justifiable if the 
objective is minimize the energy loss in the atmospheric phase of the 
maneuver. In this case one can treat energy as constant in the dynamics, and 
account for the energy loss in the performance index using the following 
integral form 


J = 


” tf C D S(l+X 2 )pV 3 /4m dt 

J o 


(17) 


Thus, E can also be regarded as an ignorable coordinate in this case. 
This will result in a reasonable approximation if the energy loss is small 
compared to the total vehicle energy. This approximation is greatly improved 
if the control solution is periodically updated to account for the present 
vehicle energy during the maneuver, which would be the case if a feedback 
(analytic) optimal control solution form was obtainable. 


2.2 Singular Perturbation Formulation 

The main approximation introduced here is that altitude and flight path 
angle dynamics can be regarded as fast compared to heading dynamics. In the 
context of singular perturbation theory, this implies a further order 
reduction to a single state model, with altitude as the control variable. To 
motivate this viewpoint, it is desirable to identify a small parameter which 
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appears as a multiplying factor on the left side of the altitude and flight 
path angle equations of motion. Currently, there is no systematic procedure 
for putting the equations of motion in this standard form. However, it is 
generally agreed that the equations of motion should always be 
non-dimensional ized as an initial step. The following transformations are 
introduced here to justify the formulation adopted in this study. 

Define the following non-dimensional variables: 



tj = (C L Sp o V c /2m<J> f )t 

(18) 

■ wv c 

Pi = P/P 0 

(19) 

■ "'"o 

r l = r/r o 

(20) 


where h is the entry altitude, p 0 is the air density at h = h , V Q - 
[v/r ]V 2 is the circular velocity at h Q , and is the final heading (final 
inclination for small changes in <f>). Also assume that for the altitudes of 
interest that r = r g . Then (9), (11-12) become: 


d^/dt^ = 

PjV^Asiny/cosy 

(21) 

edh^/dtj 

= V^siny 

(22) 


= C L SpV c V 1 Acosy/2m + V^Cl-l/V^jcosy/r,. 

(23) 


where e = C*S p o h o /2mij> f . Hence e is a small parameter for sufficiently large 
heading changes, or if 

C*Sp 1 h o /2m«i| )f ( 24 ) 

A typical calculation for a 40° plane change and h Q = 200,000 ft gives e = 
0.0043 for a vehicle with a maximum L/D of 2.3. 

The analysis in Section 3 of this report uses the original state 
variables and artificially introduces e = 1.0 as a scaling parameter. It can 



be shown that this formal procedure results in the same control solution as 
that obtained using the formulation in (21-23) in the non-dimensional ized 
variables. 

2.3 Boundary Conditions 

In [6], the sensible atmosphere is assumed to occur at h(0) = 200,000. 
ft. The starting velocity and flight path angle are derived using a deorbit 
impulse from circular orbit at h c = 100 nm, which is optimized for the 
atmospheric maneuver of interest. The initial heading angle is taken as 
zero. In the singular perturbation formulation, altitude appears as a 
control variable in the reduced problem. The optimal solution appears in the 
form 


h = h Q (E) ( 25 ) 

For comparison purposes, in this paper the starting energy is chosen to match 
that of [8] for the case of a 40° heading change. From conservation of 
energy this results in the same deorbit impulse. The initial flight path 
angle is derived from conservation of angular momentum. 

y(o) = -cos *[(r c ) (V c “A\/^)/(r s +h(o)V(o) )] (26) 

where r is the circular orbit radius and V =(v/r )^^ is the circular velo- 
C u ^ 

city. The vehicle begins the maneuver with a mass m c and, as a result of 
the deorbit impulse, the mass for the atmospheric portion is given by 

m = m c exp(-AVj/C) (27) 

The terminal condition is taken as: 

*(t f ) = > 0 (28) 

For the reentry problem there are no terminal constraints on altitude and 
flight path angle, thus their corresponding costate values are zero at the 
final time. For the noncoplanar orbit transfer problem, the final altitude 
is constrained to ensure exit from the atmosphere. 



SECTION 3 

HEADING CHANGE WITH MINIMUM ENERGY LOSS 

The objective is to minimize the energy lost in maneuvering to a speci- 
fied heading. Regarding energy as a slow variable, and altitude and flight 
path angle as fast variables, the following three state model was adopted for 
singular perturbation analysis: 



* 

di|>/dt = C L pSVAsiny/2mcosy 

(29) 

w 

edh/dt = Vsiny 

(30) 


sy/dt = C*pSV(Acosp+Mcosy)/2m 

(31) 

v# 

The objective is to minimize 



tf 



J = - I £dt 

J n 

(32) 


o 


where 

£ = -C*(l+A 2 )pSV 3 /4m (33) 

Note that in the above formulation E is approximated as constant in the 
dynamics, but that changes In energy are accounted for in the performance 
index. This approximation will later be relaxed in the subsequent analysis. 
The perturbation parameter e is introduced to signify the presence of fast 
dynamics, and is nominally equal 1.0. We seek a reduced and zero order boun- 
^ dary layer solution about e = 0, in accordance with the procedures detailed 

in [2-4]. Regarding h and y as fast states is characteristic of energy state 
analysis for fighter and transport aircraft. Therefore, we adopt the same 
framework in this analysis. Considering both h and y in the same time scale 
results in a two point boundary value problem. A feedback guidance law is 
obtained by expansion of the necessary conditions to first order [15]. In 
this regard, it should be noted that there have been some studies that have 
also considered analysis of h and y dynamics on separate time scales [16], 
which avoids linearization the boundary layer dynamics. Therefore, an 
additional set of guidance laws are possible beyond those presented here. 
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3.1 Singular Perturbation Analysis 
Reduced Problem 


Setting e = 0 in (29-31) the necessary conditions for optimality become: 


N-? 


H = A.<i> 
o </> 

o 

11 

•UJ 

1 

(34) 


V 

1 = 0 , 

Acosp = -M 

(35) 

V 


V h o * 

arg min {<!>/£} 
h,p 

(36) 




It can be shown that this results in the following reduced solution: 

X Q = (1 + ( 37 ) 

sinp o = [(1 + M*)/(l + 2M^)] 1/2 (38) 

h = arg min (v' (1 + M ) /2 }| E = const (39) 

h 

where M is the value of M for h = h . The quadrant for the bank angle in 
o ° 

(38) is resolved based on the following inequalities: 

0 < p < n/2 for M < 0 (40) 

o o 

-rr/2 < p < it for M > 0 (41) 

0 0 

It can be seen from the above solution that M plays a crucial role in the 
solution process. 

Since most of the energy is kinetic, V is weakly dependent on h for con- 
stant E. This can readily be seen from (16) where changes in h give rise to 
small changes in r. Thus, the minimization in (39) results in a value for M 
very close to zero. The interpretation is that the maneuver should be per- 
formed at an altitude where gravitational and centrifugal forces nearly 


V 
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w 


V 




cancel one another. For M small, it can be seen from (37,38) that the maneu- 
ver is performed at near maximum L/D and at near 90° of bank angle. These 
results are in good agreement with the results in [6] for the AOTV problem. 
It will also been shown in Section 3.2 that the reduced solution altitude 
profile, h (E), closely resembles the altitude profile of [6] for the case of 
large changes in inclination angle. 

Boundary Layer Problem 

Introducing the transformation x = t/z and again setting e = 0, the nec- 
essary conditions in the boundary layer are: 

h BL = * x h Vs1 " Y + v ' t = 0 

3H b( _/3Lj = 0, 3H bl / 3L 2 = 0 
where A 0 is determined in the reduced solution from (34) 

■ E°/*° * - c d v o< 1+m o) 1/2/c l 

using the solutions for A Q , y Q and h Q . In (43), and Lg represent the 
horizontal and vertical components of lift coefficient 

L« = Asiny l _2 = Acosy (45) 

which are now used as control variables in place of A and y. 

The first condition in (43) results in 

L * = (V q /V) 2 ( 1 + M 0 ) /z /cosy (46) 

where M , V are the values of M and V corresponding to h = h Q for the 
current value of E. This solution approaches the corresponding reduced 
solution as h approaches h . 

The second condition in (43) yields 


(42) 

(43) 


- .10 - 



( 47 ) 


* * * 2 
L 2 = -<C L /C D V U Y 

which can also be shown to approach the reduced solution as h approaches h Q , 
where 


>o _ 


•jf 7 it 

= C D V o W 




(48) 


The reduced solution for X in (48) follows immediately from (31) and (47) 

with 1*2 = -M q (ey = 0, in the reduced solution). The second condition in 

(48) is a consequence of the fact that h Q results from an unconstrained 

minimization of H°. Note that a£ satisfies the terminal boundary condition, 

but A 0 does not. This point will be addressed in the numerical results 
Y 

section. 

Unfortunately, evaluation of A^ needed in (47) requires the solution of 
a two-point boundary value problem. When close to the reduced solution it 
may be possible to use (48), which results in the following expression for 
flight path angle rate 

dy/dt = C^pSV(Mcosy " V q M 0 /V )/2m (49) 


For y near zero and h near h , (49) simplifies to 

dy/dt = C*pSV o (M-M Q )/2m ( 5 °) 

Use of (46) and (47) with A^ = A° results in a guidance law in feedback form, 
which we denote as the "SP1" Solution. 


Expansion of the Boundary Layer Problem [15] 

A second feedback solution can be obtained by considering an expansion 
of the boundary layer necessary conditions of (42,43) together with the state 
and costate dynamics expressed in the stretched time scale x = t/e: 

dh/dx=V siny, dA^/dx=-3Hg^/3h (51) 
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(52) 


dY/dx=C* p SV(L 2 +McosY)/2m, dA Y /dx=-9H BL /9Y 

Substituting for 1^ and L 2 from (46,47), equations (51,52) are expanded about 
the reduced solutions equilibrium conditions: 


h = h 0 (E), 


A h - 0 


= 0 


A = A° 
Y Y 


(53) 

(54) 


This results in the following linear perturbation equations: 


6h 1 


Y 1 


*'h 


6 A * 



0 

0 

K, 


0 

0 

0 

-V 


0 

K, 

C 

-K, 



Sh ' 


Y 


A 


5A 

yJ 


(55) 


where 


K : = [V 2 r 2 -(pV 2 +y)r+23p]/pVr 3 


■2p 2 /V 3 r 4 - C*SpMp/mVr 2 


K 2 = _C L Sp/2mC D V 


(56) 

(57) 


K 3 = -9 2 H BL /9h 2 < 0 


(58) 


K 4 = C*SpV 3 (l+2M 2 )/2m 


(59) 


The expressions in (56-59) are evaluated at h = h Q , and p is the scale height 
in an exponential atmospheric model. The term in (58) is complicated to 
express analytically, but can easily be evaluated numerically taking into 
account the fact that both H gL and 3Hg L /9h evaluated at h = h Q is zero. 

The eigenvalues of (55) are arranged symmetrically about the imaginary 
axis, and occur in complex conjugate pairs. In order to suppress the 



instability in the boundary layer response, the state vector in (55) is 
expressed as 

x = k^a + k (60) 

where x T = [6h, y, V and t,t are the real and imaginary parts of the 

eigenvectors associated with the stable eigenvalues. Knowing 6h and y, it is 
a simple exercise to solve for kp kp A^ and <5Ay. Then in (47) can be 
evaluated for 

A = A + 5A ( 61 ) 

y Y Y 

Equations (46,47) and (61) thus constitute a second feedback guidance law, 
which we term the "SP2" Solution. 


Modeling Energy Rate Dynamics 


If energy rate is modeled in the dynamics, the reduced model becomes a 
two-state problem, and the performance index is modified to minimize -E(t^). 
The Hamiltonian in this case is 

H o = A J» + A e £ = 0, A E (t f ) = -1 (62) 

This gives rise to a two-point boundary value problem in the reduced solu- 
tion. However, an approximate integration of A £ is possible in this case, 
based on the known properties of the optimal solution. Using (62), it is 
easy to demonstrate that 

dA £ /dt = -9 H q /9E * A e C d Spi/(1+A )/2m (63) 

Thus, 

dA E /dip = (2C*(l+A 2 )cosy/C*Asinu)A E = (2C*/C*)A E (64) 



where the approximations A = 1, p - tt/ 2 , y = 0 have been employed. Integra- 
tion of (64) results in 

A e (i|>) = -exp{-2Cp(ii> f -<|>)/C L } ( 65 ) 


Comparing the Hamiltonian expressions in (34) and (62), it can be seen 
that modeling E as constant in the dynamics amounts to the approximation A £ = 
-1. Equation (65) represents an improvement, but the approximation A g = -1 
is apparently accurate for high L/D vehicles. So far as its effect on the 
reduced solution, (37-41) remain the same. The reduced costate solutions 
become: 


A° = 


* ? ? 

-Wo (1+I V > 


1/2 


/C. 


A° = 0, 


A 0 = A r c n v ^ Mo/ 0 !" 

Y E D 0 0 L 


( 66 ) 

(67) 


Note that the $ and y costate solutions are now simply multiplied by A g . The 
boundary layer solution for l* in (46) remains the same, but (47) becomes 

L 2 = W C D v2) V X E (68) 

Thus the SP1 control solution, which uses A = A° remains unchanged when E is 
modeled in the dynamics, since A g is canceled when A^ from (67) is substi- 
tuted in (68). The SP2 solution, on the other hand, is affected in that 
several of the matrix elements in (55) are changed. In particular, K 2 and K 4 
are divided by A £ and H gL used in the computation K g becomes 


H 8L = V + V s1nY * V + <69> 

We will refer to the control solution obtained with these modifications as 
the "SP3" Solution. 


3.2 Numerical Results 


A numerical study was performed to evaluate the performance of the three 



guidance laws derived in the preceding section. The parameter values, chosen 
to approximate the vehicle studied in [5], are as follows: 

C Q = .032, S = 125.8 ft 2 , K = 1.4, m = 331.5 slugs 
0 

The initial conditions were chosen to match the 40° heading change maneuver 
of [5], where the sensible atmosphere was defined to begin at an altitude of 
200,000 ft. The corresponding entry velocity and flight path angles are V(0) 

= 25,945 ft/s and y(0) = -.148°. A simple exponential atmospheric model was 
defined using the standard atmospheric data for air density at altitudes of 
10 5 and 2x10^ feet. All of the comparisons are for a 40° heading change. 

We first illustrate the validity of the singular perturbation 
formulation by comparison of the reduced solution altitude profile obtained 
from the use of (39), with the optimal solution for the 40° plane change 
problem. The optimal solution was calculated using a multiple shooting 
method described in [17]. Figure 1 illustrates this comparison. The optimal 
solution altitude profile satisfies a terminal constraint that h(t^) - h(0), 
needed for a typical A0TV orbit plane change. It can be seen that the 
reduced solution altitude profile closely follows the optimal altitude 
profile, with the exception of satisfying the initial and final values of 
altitude, which are lost in the reduced problem formulation (altitude is a 
control like variable). Clearly, the reduced solution can be used for the 
initial phase of an A0TV plane change maneuver, but a large terminal boundary 
layer correction is needed for the exit phase. This aspect will be addressed 
in Section 4 of this report. 

We next consider the performance of the guided solutions for the reentry 
problem. Figure 2 compares the reduced solution altitude profile obtained 
from (39) with the SP1 and SP3 guided solution profiles, which on this scale 
are identical. A similar comparison is given for the SP2 guided solution in 
Fig. 3. Note that the reduced solution provides a reasonable altitude 
profile except at high energies near the initial time. However, this region 
of the solution is of little interest since the air density is negligibly 
small. In any case, it is not physically possible to follow this profile 
(recall that h is used as a control variable in the reduced solution). Thus, 
it was decided to maintain A = 1 and Acosy = -M (y = 0) until h Q (E) falls below 
the current altitude. 



In order to evaluate the optimality of these solutions, an optimal solu- 
tion was numerically computed using the four-state model in (8-11) to define 
the dynamics. The SP1 guided solution was used as an initial guess for the 
state time histories, and the reduced solutions in (44), (54), and (65) were 
used as an initial guess for the costate time histories. The solution 
converged to a relative precision of 10 in eight iterations. The value of 
the Hamiltonian was constant and essentially zero considering the relative 
precision accuracy required for convergence. This served as an independent 
check on the accuracy of the solution. 

Figure 4 compares the optimal altitude profile with the profiles that 

result from the three guided solutions. Note that the optimal solution dips 

slightly more into the atmosphere near the end, and consequently results in 

slightly decreased flight time. The corresponding control time histories and 

heading profiles are compared in Figs. 5-7. Note that in Fig. 5, the optimal 

bank angle at the final time is 90°, which follows from (47) and the fact 

that A (t £ ) = 0. In the context of singular perturbation theory, this gives 

Y T , 

rise to a terminal boundary layer which must be solved backwards in time. 

Since this was ignored in our analysis, the guided solutions approach the 

condition in (48) instead. This explains the departure in the altitude 

profiles of Fig. 4. It is apparent that this effect is a minor one. In any 

case, the dip in the optimal profile may not be desirable from a practical 

standpoint. 

It may be somewhat surprising at this point that the SP1 and SP3 solu- 
tions are nearly identical. However, recall that the SP1 solution is not 
sensitive to the approximation that E is modeled as constant in the dynamics. 
The SP3 solution corrects the SP2 solution for this modeling approximation, 
and results in essentially the same solution as the SP1 solution. This fact 
also justifies the use of approximation A^ = A^° in the SP1 solution. Table 
1 compares the energy loss for all the solutions, and shows the three guided 
solutions produce essentially optimal performance. The energy loss for the 
SPl and SP3 solutions is indistinguishable from the optimal solution to four 
significant places, while the energy loss for the SP2 solution is .06% 
greater. 



SECTION 4 

AEROASSISTED ORBIT TRANSFER 


The preceeding section has demonstrated that singular perturbation 
theory can be used as an effective tool in model order reduction for the 
problem of hypervelocity heading change within the atmosphere with minimal 
energy loss. The extension to aeroassisted orbit transfer with inclination 
change requires that a contraint be imposed on the terminal altitude to 
insure that the vehicle exits the atmosphere when the heading change (which 
approximates the inclination change) is achieved. This introduces a 
difficult terminal boundary layer problem, which to date we have not been 
able to solve in closed form, since it requires an analytic integration of 
the boundary layer equations. In [5-9] this same problem arises, but in a 
different context. In these studies, the states and co-states associated 
with a four-state model are analytically integrated with the assumption that 
Loh's term, M(h,V), is constant or piecewise constant over the trajectory. 
While this is true for a large portion of the maneuver, M undergoes a large 
variation near the end. Consequently, the guidance algorithms resulting from 
these studies require large angles of attack near the end of the maneuver to 
compensate for this variation. This drawback was subsequently removed in 
[10], however, the resulting guidance algorithm requires that a complex 
quadrature be performed at each guidance interval. Moreover, this approach 
is not readily extended to the case where heating rate is constrained. 

The essential problem in all these approaches lies in the fact that air 
density decreases exponentially as altitude increases, and corrections to 
satisfy terminal constraints must be accurately predicted while the vehicle 
is at lower altitudes. 

In this section, the simplest guidance algorithm presented in Section 3 
is used for the initial portion of the maneuver, and a predictor/corrector 
type algorithm is presented for the terminal maneuver. The 
predictor/corrector algorithm relies" on bank angle control alone, and thus 
avoids the problem of large angle of attack. The form of the guidance 
algorithm was chosen to closely approximate the known properties of the 
optimal solution, while permitting an accurate integration of the equations 
of motion. The predictor/corrector algorithm provides the information needed 



to both initiate the terminal maneuver, and to guide the vehicle in closed 
loop fashion. The availability of a closed loop guidance algorithm will be 
of paramount importance in future studies addressing the effect of 
atmospheric anomalies, and of course, for real time implementation. 
Comparisons are made to numerically optimized trajectories for a range of 
orbit plane angles to demonstrate the near optimality of the complete 
guidance algorithm. 

4.1 Guidance During the Reentry Phase 

Guidance during the reentry phase is based on the zero order reduced and 
boundary layer solution referred to as the SP1 solution in Section 3. This 
consists of first calculating the reduced solution in (39), and then 
calculating the horizontal and vertical components of the normalized lift 
vector using (46,47), with = X^° as given by (48). 

4.2 Guidance During the Exit Phase 


The reentry phase guidance algorithm does not satisfy the terminal 

constraint on h, which was lost in the reduced formulation. The terminal 

boundary layer necessary conditions are identical in form to those for the 

initial boundary layer. However, the solution asymptotically approaches the 

reduced solution backwards in time, starting from the terminal constraint on 

altitude. In addition, the change in * during the terminal maneuver must be 

accounted for to insure that both terminal constraints are met 

simultaneously. This requires an anlytical integration of the state and 

costates. To circumvent this problem, a predictor/corrector 

guidance law was developed based on the known properties of the optimal 
* * 

solution: A =1,0, y 

To simplify the problem of integrating the forward motion, it is assumed 
that the nominal exit maneuver consists of two regions. In the first, a 
constant (negative) bank angle perturbation is used to increase the flight 
path angle. This is followed by a second region in which y is influenced 
only by gravity. A bank angle correction is computed throughout the maneuver 



based on the predicted heading error at the final altitude. Constraints are 
enforced to insure continuity at the junction of the two regions, and 
satisfaction of the terminal constraint on altitude. The maneuver is 
depicted in Fig. 8. 

During terminal guidance, we maintain A = 1.0, and modulate the bank 
angle according to the following equations: 


y = y Q - 5y Q + 5y = y + <5y 


where in Region 1: 


y = cos" 1 ( -M) , 6y = const. > 0 
0 0 


and in Region 2: 


y Q = 90% 6y 0 = 0 


(70) 


(71) 


(72) 


The nominal trajectory can be analytically predicted for 5y=0, and using the 
approximations: 


COS y = -M + 6y o siny Q , sin y = y 


(73) 


the details of which can be found in [18]. The predicted heading change 
(a<|> ) for the exit maneuver is calculated at each integration step along the 
trajectory, and is used to both initiate the maneuver (when = ^ go ) , and 
to track the terminal constraint using a simple proportional control law, 


6y = k((p 


go 


A V 


(74) 


Region 2 guidance is initiated when the present altitude satisfies the 
continuity constraint at the junction of the two regions. 

4.3 Numerical Results 


A numerical study was performed to evaluate the performance of the 
sub-optimal guidance algorithm described in the preceding section, using the 
same vehicle data as presented in [6] and in Section 3. Fig's 9 and 10 
illustrate the guided altitude profiles and the corresponding bank angle 
profiles obtained for heading changes up to 40°, in increments of 10°. The A 
profiles are very close to 1.0 throughout for all the maneuvers, and are not 



illustrated. These results were generated for 6 u q = 25°, k = 50, however it 
was found that the general character of the solutions did not change as these 
guidance parameters were varied. Note from Fig. 10 that, following the 
initial perturbation, bank angle continues to decrease in Region 1. This is 
due to the variation in M that takes place as the altitude departs from the 
reduced solution profile. After completing Region 1, the bank angle remains 
very close to 90°, or in other words, 5p in (24) is very close to zero, 
indicating the accuracy of the prediction algorithm. 

In order to evaluate the optimality of the guidance algorithm, optimal 
solutions were obtained numerically using the multiple shooting algorithm 
[ 17 ]. Fig's. 11 and 12 illustrate the optimal altitude and bank angle 
profiles. The most remarkable characteristic in these solutions is that the 
final time is nearly independent of the final heading. Also note that the 
final bank angle is always 90°, which is required by the necessary 
conditions. The corresponding A profiles are shown in Fig. 13, which 
verifies that the optimal solution lies close to A = 1. The flight path 
angle histories for all the cases in Fig. 9 and Fig. 12 are quite small, and 
close to zero at the final time. 

Despite the fact that the final times are considerably different, the 
guided solution performance is not far from optimal. Table 2 illustrates the 
near optimality by comparing the energy loss of the guided solutions with 
that of the optimal solutions. 



SECTION 5 

CONSTRAINED AEROGLIDE 

The constrained aeroglide problem is treated in this section. This is 
essentially the same problem addressed in Section 4 with the addition of a 
heating rate constraint, which amounts to a first order state variable 
inequality constraint. 

5.1 Problem Formulation 

For this problem it is necessary to consider the system of -equations in 
(2-6), since the approximation that 0=0 is no longer valid. In this case, 
9 remains an ignorable coordinate, and (7) is used to define the change in 
inclination angle. It was also necessary to initiate the maneuver at a much 
higher altitude, since the starting condition for the results in Section 4 
has a fairly high value of heating rate. Only numerically optimized 
solutions were considered, using the same multiple shooting method that was 
employed for the unconstrained solutions. First the touch point solution was 
considered, and the range of maximum allowable heating rate over which this 
solution applies was obtained for the same vehicle dynamics considered in 
Section 4. Then, the constrained arc solution was attempted for lower values 
of the heating rate limit. However, only a narrow range of solutions were 
found for this case. Unfortunately, all of these solutions violate the 
practical limits for heating rate. Hence, the optimal aeroglide problem 
remains an open research issue. It may be that very complex behavior results 
for lower values of heating rate limit, such as multiple touch point 
solutions combined with portions of a constrained arc. 

The expression for heating rate employed in this study was: 

dQ/dt = 17600[p/pJ 1/2 [V/V ] 3 ' 15 (BTU/sec/ft 2 ) (75) 

S 3 


5.2 Numerical Results 

The vehicle characteristics are the same as those given in Section 3.3. 
The initial and terminal condition data are as follows: 



h ( 0 ) = h(t f ) = 365,000 ft, V(0) = 25,746 ft/s, y(0) = 0.55° 

1(t f ) = 18° 

The initial h, V and y correspond to the same deorbit impulse employed in 
Sections 3 and 4. 

Figures 14-16 illustrate the h, V and i profiles for the unconstrained 
solution and for the touch point solution corresponding to a maximum heating 
rate of 600 BTU/sec/ft 2 . As the heating rate limit is decreased, the minimum 
altitude and flight time increase. Figures 17 and 18 show the corresponding 
control time histories. For the touch point solution, there is a slight jump 
in the normalized lift coefficient at the touch point. However, this is 
likely due to numerical inaccuracy in satisfying the touch point conditions. 
The bank angle history is continuous. Figures 19 and 20 illustrate the jumps 
that occur in the A^ and A ^ at the touch point. All other costates are 
continuous. All of the state and costate histories can be found in [18]. 
Figure 21 compares the heating rate profiles. The peak heating rate for the 
unconstrained case was 729.3 BTU/sec/ft 2 . In comparing these results it was 
found that the increase in energy loss due to the heating rate limit imposed 
was negligible. However lower heating rate limits will likely result in much 
greater energy loss. 

Attempts to decrease the heating rate limit below 600 BTU/sec/ft 
resulted in convergence failure. This failure was abruptly encountered at a 
heating rate limit of 599.3 BTU/sec/ft 2 . Thus we suspect that either a 

conjugate point is encountered, or the nature of the optimal solution changes 
to either a constrained arc case or a multiple touch point case. The 
constrained arc case was investigated, and after many attempts (including 
reformulation in nondimensional variables) we were not able to obtain a 
converged solution. 


- 22 - 



SECTION 6 

CONCLUSIONS AND RECOMMENDATIONS 


The application of singular perturbation methods to optimal control 
problems related to aeroassisted orbit transfer vehicles has been 
investigated in this study. Two closely related problem formulations were 
addressed: optimal reentry of a hypersonic gliding vehicle, and optimal 

orbit plane change with aeroglide. In addition, an attempt was made to 
obtain numerical solutions for optimal aeroglide orbit plane change subject 
to a maximum heating rate constraint. 

6.1 Conclusions 

The major conclusions resulting from this research effort are as 
follows: 

(1) Singular perturbation theory using energy state approximations can 

be used to reduce the model order to a single state equation, and a 
closed form solution for the reduced problem can be readily 
obtained. The solution for the reduced problem reasonably 

approximates the full order optimal solution, except near the 
initial and final conditions. By non-dimensional izing the 
equations of motion a singular perturbation expansion parameter can 
be identified that depends on vehicle parameters and the required 
heading change. 

(2) For the reentry problem, three guidance algorithms were derived, 

all of which are nearly optimal in terms of minimizing the energy 

loss for the maneuver. No terminal boundary layer correction was 
required for this problem formulation. 

(3) A constraint on terminal altitude, required for the noncoplanar 

orbit transfer problem, results in a difficult boundary layer 
problem for which we were not able to obtain a tractable analytic 
solution. However, the optimal terminal maneuver was approximated 



using a sub-optimal predictor/corrector guidance law. In general, 
optimization of the terminal maneuver (to satisfy terminal 
constraints) is not as critical as optimizing the initial (reentry) 
portion of the maneuver. 

(4) The problem of optimal orbit plane change subject to a heating rate 
constraint results in a touch point extremal solution for a high 
(but narrow) range of heating rate limit. These trajectories and 
the corresponding control histories are similar to the 
unconstrained solution, with negligible increase in energy loss. 
This will not be true for lower values of the heating rate limit, 
where the constrained trajectory may ride the constraint boundary 
and/or contain multiple touch points. Unfortunately, the multiple 
shooting method could not be successfully used to find extremal 
solutions corresponding to the first order necessary conditions 
associated with this problem. 

6.2 Recomnendations 

Based on the results of this research effort, the recommendations for 
further work along this line are: 

(1) The numerical study of optimal aeroassisted orbit transfer with 
aerogl ide, subject to a heating rate constraint, should be 
completed over a more practical range of heating rate limit. Along 
this line, alternative formulations of the necessary conditions, 
such as the transformation method in [19], should be explored. 
Another approach is to use a parameter optimization based method 
such as that employed in OTIS [20]. In this case, the constraint 
would be enforced indirectly through a penalty function approach. 
The resulting profiles are of interest for comparison to aerocruise 
solutions [12], and for further developing near optimal guidance 
algorithms based on singular perturbation theory. 

(2) Further research is needed to develop methods of solving the 
terminal boundary layer problem associated with terminal 



constraints. One such alternative, investigated in this study [18] 
(but not reported here), is a formulation in which the terminal 

altitude constraint was satisfied as a part of the reduced 

solution. This avoids the terminal boundary layer problem 

altogether. 

(3) Parametric studies of the effect of initial and final orbit 
altitudes on the optimal orbit plane change maneuver should be 
conducted. This includes the possibility of optimizing the deorbit 
impulse, and including the use of multiple impulsive maneuvers for 
performing a part of the orbit inclination change outside the 
atmosphere. 

(4) Robustness of the guidance algorithms resulting from this study was 

not investigated. In particular, the effect of uncertain 

atmospheric conditions at high altitudes should be evaluated. 

(5) Extensions to problem formulations suitable for aerocapture and 
orbit plane change for future Mars missions should be explored. 

6.3 Publications 

The journal and conference publications that resulted from this research 
effort can be found in [18] and in [22-26]. 
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Figure 1. Comparison of the reduced solution altitude profile 
with an optimal profile for a 40° heading change. 








Figure 2. Comparison of the SP1 and SP3 guided altitude profiles 
with the reduced solution. 
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Figure 4. Comparison of the guided altitude profiles 
with the true optimal solution. 
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Figure 7. Heading profiles. 




Figure 8. Depiction of the exit phase maneuver. 
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Figure 9. Guided solution altitude profiles. 








Figure 10. Guided solution bank angle profiles 
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Figure 11. Optimal solution altitude profiles. 






Figure 12. Optimal solution bank angle profiles. 
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Figure 13. Optimal solution normalized lift coefficient profiles. 
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Figure 14. Altitude profiles for the unconstrained and touch point 
solutions. 
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Figure 15. Velocity profiles for the unconstrained and touch point 
solutions. 
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Figure 16 . Inclination angle profiles for the unconstrained and touch 
point solutions. 
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Figure 17. Normalized lift coefficient profiles for the unconstrained 
and touch point solution. 
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Figure 20. Velocity costate profiles for the unconstrained and touch 
point solutions. 



Table 1 

Comparison of Final Energies for the Reentry Problem 


Guidance 

Final Time 

Efxl08(ft 2 /sec2) 

AE 

Optimal 

358.6 

-4.813 

1.510 

SP1 

397.0 

-4.813 

1.510 

SP2 

415.8 

-4.814 

1.511 

SP3 

398.0 

-4.813 

1.510 


Table 2 

Comparison of Energy Loss for the AOTV Problem 


Heading 

Change 

Optimal 

Solution 

xlO 7 

Guided 

Solution 

xlO 7 

%Error 

10° 

4.713 

4.971 

5.47 

20° 

8.858 

9.109 

2.83 

30° 

12.38 

12.63 

2.04 

40° 

15.40 

15.66 

1.69 
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